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Abstract. We derive the steady state properties of a general directed "sandpile" 
model in one dimension. Using a central limit theorem for dependent random 
variables we find the precise conditions for the model to belong to the universality 
class of the Totally Asymmetric Oslo model, thereby identifying a large universality 
class of directed sandpiles. We map the avalanche size to the area under a 
Brownian curve with an absorbing boundary at the origin, motivating us to 
solve this Brownian curve problem. Thus, we are able to determine the moment 
generating function for the avalanche-size probability in this universality class, 
explicitly calculating amplitudes of the leading order terms. 
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1. Introduction 

Sandpile models have played an important role in developing our understanding 
of self-organized criticality [HI IH El El E|. One important notion is that of 
universality, the idea that quantities such as critical exponents and scaling functions 
are independent of microscopic details of the model. This has been studied in 
the context of individual models, but few have determined general conditions for 
models to belong to a particular universality class [U El E] • In the following, we 
present details of the solution of a general directed one-dimensional sandpile model 
introduced in [TU] which is a generalisation of a model studied in (TJ EJ- We use a 
central limit theorem for dependent random variables ^T] to determine the precise 
microscopic conditions for scaling of the moments of the avalanche-size probability. 
We also argue that there is an n-dependent crossover length £ n , such that for systems 
with size L <C £ n branching process behaviour is observed. 

The avalanche size statistics are calculated by mapping the model to the 
problem of finding the area under a Brownian curve with an absorbing boundary at 
the origin, that is, if x(t) is the trajectory of a Brownian curve such that if x(t') = 
for some t' > then x(t > t') = 0. In the large L limit, the avalanche size statistics 
are identical to those for the area under the Brownian curve after a "time" equal to 
L; A = x(t)dt. This motivated us to calculate the moment generating function 
for this area, which is an interesting problem in its own right as there have been some 
recent interest in physical applications of the statistics of the area under Brownian 
curves H2 E31 EE3 ■ 

2. Definition of The Model 

The model we study is on a one-dimensional lattice of length L where each lattice 
site, i = 1,2, ... ,L, may be in one of n states. The state of site i is denoted z i: 
which may take values 0, . . . , n — 1 and this is interpreted as the number of particles 
on site i. 

At the beginning of each time step a particle is added to site 1: z\ — > z\ + 1. 
This site may topple a number of times, each toppling redistributing one particle 
from site 1 to site 2: Z\ — > z\ — 1, z^ — > + 1. When site 2 receives a particle it may 
undergo topplings, redistributing particles to site 3, which in turn may topple, and 
so on until either a site does not topple, or site L topples where the redistributed 
particles leave the system and the time step ends. The avalanche size, s is the total 
number of topplings which occur during a single time step. The toppling rules are 
therefore defined through choosing the probability that a site with z particles will 
topple so many times upon receiving a particle. 

The only restrictions on the topplings are: (i) Zj must remain in the range 
[0, n — 1]. For instance, a site with Zi = 2 may not topple more that three times 
when receiving a particle. Moreover, a site with z% = n — 1 which receives a particle 
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must topple at least once., (ii) When site i topples it redistributes exactly one 
particle to site i + 1 only, (iii) The toppling rule is homogeneous and obeys a 
Markov property in that the probability that site i topples s$ times depends only 
on Zj. In fact, the requirement on the homogeneity can be relaxed and all the 
following is trivially extended to inhomogeneous toppling rules. Each site will then 
have a different stationary state, but provided that the remaining constraints are 
obeyed, the scaling of the avalanches will remain unaltered, (iv) There must be some 
probabilistic element to the toppling rules. To be precise, there must exist at least 
one value of z such that the number of topplings a site in this state undergoes is 
non-deterministic. This last restriction discounts purely deterministic toppling rules 
which lead to trivial dynamics. 

Self-organised criticality (SOC) is associated with a stationary state where the 
avalanche-size probability, P(s; L), which is the probability of observing an avalanche 
of size s in a system of size L, obeys simple finite-size scaling, 



P(s; L) = as- TC S (s/bL A ) for s > 1, L > 1 



where a and b are non-universal constants, r and A are universal exponents and £f 
is a universal scaling function. It can be shown that if lirn^oo Jy hL A u k ~ T< 3(u)du 



exists and is non-zero for k = 1,2, 



then 



s=l 

oc L^ k+1 - T \ 



(2) 

Hence, for L ^> 1, the scaling of the moments, (s k ) l with L is universal and we can 
determine the universality class of a model by calculating the exponents A and r 
from the scaling of these moments. 

To formulate the model, we use a Markov matrix representation, which is an 
extension to the work in [2|. The configuration of a system is {z\, z 2 , ■ ■ ■ , zl} for 
which we shall use the shorthand notation {z{\. In order to construct a Markov 
matrix representation, we first consider a representation for the configuration of a 
single site, that is, a system of size L = 1. Such a system can be in one of n stable 
configurations, z — 0, 1, . . . , n — 1, and therefore we require an n- dimensional vector 
space to represent a probability measure over these configurations. We therefore 
construct n-dimensional left and right basis vectors, {(e 2 |} and {|e z )} where 



z times 



(0,...,0, 1,0, ...,0) and \e z ) 



0-f 

1 


v o y 



z times 



(3) 
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such that (e z \e z >) = 5 Z)Z t. The configuration of a single site with z particles is then 
represented by the basis vector \e z )i = \e z ), where the subscript 1 reminds us that 
it is the configuration of 1 single site. The n 1 '- dimensional vectors representing 
configurations [zi] for systems with L sites, \e{ Zi })i, are constructed from these 
basis vectors 

\e {zi} ) L = \e Zl )i ® \e Z2 )! <g> . . . <g> \e ZL )i (4) 

where ® is the usual tensor product. 

A state vector \Pt)L is the weighted sum of basis vectors, 

\Pt)L = ^w\ Zi} \e {zi} ) L , (5) 
{*} 

where the weights ^L.} are the probabilities that the system is in the configuration 
{z^ at time t, that is w\ Zi \ = ( e {zi}\Pt)- Note that the normalisation condition 
requires 



{Zi} {Zi} 



ir;. , = $> {jN} |P t >£ = l- (6) 



The system is evolved by applying operators to the state vector \Pt)i- We 
define a toppling operator Gl which adds a particle to site 1 of a system of size L 
and carries out all the topplings: 

\Pw)L = G L \P t ) L , (7) 

making it clear that \Pt)i is a Markov chain. Using the |e z ) representation the 
toppling operator is an n L xn L matrix where {e^^Qi^e^) i is the probability that 
adding a particle to a system in the configuration \e^ Zi y) i results in the configuration 
\z{z'.})l after topplings. See Appendix |Appendix A| for an explicit representation of 
Q>l for a system with n = 2. The steady state, \0)l is defined as the state which 
is invariant under application of the toppling operator G^. It is therefore the right 
eigenvector of G^ with eigenvalue 1, 

Gx|0>x = |0>l. (8) 
The corresponding left eigenvector, (0\ L satisfies 

(0| L G L = (0\ L . (9) 

Since particle number is conserved in the bulk, the sum of each column in Gl equals 
1. This leads to 

(1,1,...,1)G L = (1,1,...,1), (10) 

identifying the left eigenvector of the toppling operator (0|l = (1,1,...,1). 

In order to calculate the moment generating function for the avalanche- size 
probability it is convenient to define the operator Gl{x) such that 

QL,m(x) = (0\ L [G L (x)} m \0) L (11) 
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is the avalanche size moment generating function over m time steps and 

<«*>l= =Qg) (12) 

x=l 

gives the fcth moment of the avalanche-size probability in a system of size L. The 
toppling operator is then 

G L = G L (1). (13) 

To illustrate how to construct G L (x), consider the action of adding a particle to site 
i = l. First, we look at the number of particles on that site, zi, and determine 
the number of times the site topples using the probabilities given by the toppling 
rule. Hence, we require n + 1 matrices of dimension n, denoted which will act on 
the first site and remove k particles, multiplying the final state by the probability 
that this toppling took place. We then need to redistribute these particles to site 
% = 2, which we achieve by acting (Gi_i) fc on the remaining L — 1 sites. This works 
because adding to site 2 of a system of size L is equivalent to adding to site 1 of a 
system of size L — 1 (note Go(x) = 1). Finally, we multiply the remaining state by a 
factor x k , which marks the state as having toppled k times, which gives moments of 
the avalanche size upon differentiation. This leads us to write the general toppling 
operator 

n 

G L (x) = 5> fc [1 <g> G L ^(x)] k S k <g> (14) 

fc=0 

where 1 is an n x n identity matrix and A® N = A <g> A <g> . . . ® A, N times. The 
restrictions, (i)-(iv), on the model give [Sjjjj > for j = i + k — 1 and equal to zero 
otherwise, that is, 

Sfc = ^ \e z +k-i)iS Z) z+k-i(e z \i (15) 

z 

where the sum is over all z which satisfy both < z < n — 1 and < z + k — 1 < n — 1 
and S z>z+ k-i is the probability that a site with z particles topples k times on receiving 
a particle. Note that particle conservation requires Y^jZl &zj — 1 

3. Stationary Properties 

In this section we find the steady state, |0)l, which is the eigenvector of Gl with 
eigenvalue 1. Consider the single site operator, 

n 

G 1 (x) = J2* k Sk. (16) 

fc=0 

We begin by finding the eigenvectors and eigenvalues defined by 

(A i (x)|G 1 (a;) = A i (A i (a;)| (17a) 
G 1 (a;)|A i (a;)) = A i |A i (a;)) (176) 
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where % takes values from to n — 1. From the properties of the and the 
normalisation condition we find 

1 



(Ao(z)l = 



1 1 

1 

x 



satisfies 



{X {x)\^x k S k = x (X (x)\ 



(18) 



(19) 



k=0 



and so is the left eigenvector of Gi(x) with eigenvalue Ao = x. The corresponding 
right eigenvector therefore satisfies 



Gi(x)|Ao(a;)) = ^VS fe |A (:r)) = x\X (x)), 



(20) 



k=0 



which determines the precise form of the eigenvector. As the eigenvectors must be 
normalised, (X (x)\Xq(x)) = 1, we may write 



|A (aO) 



/ Pox 71 ' 1 \ 

PlX n ~ 2 



(21) 



V Pn~\X° ) 

where p z is the probability that a site contains z particles in the stationary state 
and Y1™ZqPz = 1- We cannot, however, determine Pi any more precisely without 
details of the and these will have to be calculated separately in each case. 

If the matrix Gi is a regular Markov matrix, that is, there exists an integer 
N > 1 such that [Gf]jj > for all then we have found the unique stationary 
state of the single site operator, 

/ Po \ 

Pi 



|0> = |A (1)> = 



(22) 



\ Pn-1 j 

In the following we shall always assume that Gi is regular. The discussion of 
the necessary and sufficient conditions for a regular Gi is non-trivial and we shall 
not discuss it in detail. We shall simply note that this requirement, along with 
restrictions (i)-(iv), still leaves an abundance of choice for the toppling rules. For 
instance, it is easy to demonstrate that any tridiagonal matrix with positive definite 
elements is regular. Hence, any toppling rule which always allows a site to topple 
zero times, once or twice on receiving a particle (with the usual exceptions for z = 
and z = n — 1) will automatically lead to an acceptable toppling rule. 

Now that we have found the stationary state of the single site operator Gi, 
we shall proceed to determine the stationary properties of the full operator Gl by 
induction. We introduce the notation \\i(x)) L as the ith eigenvector of Gl(x) such 
that 

G L (x)\X i (x)) L = X Lti (x)\X^ L , (23) 
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where \l,i(x) is the ith eigenvalue of Gl(x). We now make the ansatz that these 
eigenvectors may be expressed 

|A i+ j(a:)>L = |Aj(a:AL-i,i(a;))>i <g> |A^(x)> L -i. (24) 
Operating on the left hand side with Gl(x) we find 

n 

G L (x)\X l+j (x)) L = J2( xX L-^) kS k\K+j(x)} L (25) 

k=0 

= Gi(xAL-i,i)|Aj(a:Ai ( _i >i (a:)))i <g> \\(x)) L -i (26) 
= x\ j (x\ L - 1 J\\ i+j (x)) L (27) 

So, we find that \\ i+ j(x))L is indeed an eigenvector of Gl(x) with eigenvalue 
^L,i+j = ^j(x\L-i,i)- Hence, by induction, and recalling that we have assumed 
Gl is regular, the unique steady state is 

\0)l = \0)i ® |0> L _x = |0)f L . (28) 

This is a product state, which means that it has no spatial correlations and indeed we 
find that avalanches in small systems are uncorrelated, leading to branching process 
behaviour. However, we will show that in larger systems temporal correlations 
develop which bring the avalanche behaviour away from the branching process to 
that characterised by the area under a Brownian curve. 



4. Toppling probability distribution 

In this section we define the toppling probability distribution and determine some 
of its properties which will be used in the next section. The toppling probability 
distribution, P(s;L,m) is defined as the probability that a system of size L in 
the stationary state undergoes a total of s topplings on receiving m particles. In 
principle, it may be calculated from the moment generating function 



P(s;L,m) = ^(0\ L [G L (x)r\0) L 



(29) 

although this is rarely a simple task in practice. Recall that we only consider toppling 
rules obeying the restrictions (i)-(iv) with a unique stationary state. First, we show 
that P(s;L,m) has a mean value 

oo 

QH = J2 sP(s;L,m)=mL, (30) 

s=0 

which is what we would expect by considering conservation in the stationary state 
since every particle that enters the system must leave through the open boundary. 
We write down the equation for the first moment 

Qil=^(0\L[G L (x)] m \0) L ^ 

n 

= mJ2(0\ikS k \0) 1 (l + Q ( ll 1<1 ). (31) 

k=0 
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Multiplying (J20j) on the left by (Ao(x)| and differentiating, we find 

(\ (x)\Y,kx k - 1 S k \\ (x)) = l (32) 



k=0 



and so (0|i Y2=o fcS fcl°)i = 1 giving Q% m = m(l + QlI^), with Qft = 1, which 
has the solution Ql m = mL. 

Next, we show that the avalanche probability may be factorized. We define 
P(s, t; m, 1, L) as the joint probability that a system of size 1+L, which has received 
m particles, undergoes s topplings in the first site and t in the remaining L sites. 

1 d s 1 d 1 ,„, 

(33) 



P(s, t; m, hL)=^ ^^^r(°U+i [Gl+i(xi, x L )] m |0) L+1 



x'i=0x'l=0 



where 

n 

G L +l(Xl, X L ) EE X l S k ® • (34) 

fc=0 

By expanding the bracket in (}3*3*j) and carrying out the differentiation with respect 
to xi, we find 

P{s,t;m,l,L) = (0\iJ2 S [Y t k * ~ s J fl S ^l )i 

{ki} \ i J i=0 

1 d* 



t! dx l L 



{0\ L (G L (x L )) s \0} L (35) 

x=0 

where S(x) is the Kronecker delta, 6(0) = 1 and <5(x) = for 0. Identifying the 
first scalar product as simply the probability that a site charged m times topples s 
times, P(s; m, 1), we have 

P(s, t- m, 1, L) = P{s; m, l)P(t; s, L). (36) 

This is simply a statement about the fact that the directed nature means that sites 
% = 2 . . . L do not have any influence on site i = 1. Hence, if we consider Sj, which 
is the number of times site i topples during a particular avalanche, this result tells 
us that in the stationary state the sequence s±, s 2 , s 3 , . . . , Sl for a single avalanche 
forms a Markov chain. This will become important later on when we come to map 
the avalanche size s = Y2i=i s i t° the area under a random walker. 

We now derive three important results for the single site toppling probability 
distribution, P(s; 1, m), which we will use later to make the mapping of avalanches 
to the area under a Brownian curve more rigorous. First, we show that the range of 
s for which P(s; 1, m) has support has an upper bound equal to 2n — 1. Second, we 
show that P(s; l,m) has a stationary distribution for large m, 

lim P(s; l,m) = VVp (m _ s)+2 (37) 

m. — >oo < ■* 
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where the sum is over all < z < n — 1 satisfying < m — s + z < n— 1. Note that 
this is a function of m — s only. Finally, we consider the width of P(s; 1, m) around 
its mean value, 



-= yy g ,,,,2 

and show that Q\ ' m > for all m. These results lead to a further result that the 
width is finite and non-zero for all m, and approaches a constant for m — > oo, which, 
again, will become important later when we map the problem to a random walker. 

The first of these results follows immediately from bulk conservation of particles: 
Consider a site with z particles, which is charged m times. After s topplings have 
taken place it will have z' = z + m — s particles. Since both z and z' must lie between 
and n—1, P(s; 1, m) may only have non-zero values for m — n + 1 < s < m + n — 1. 
Hence, P(s; 1, m) = for \s — m\ > n — 1 since such topplings are always impossible 
and Q^ln < (n — l) 2 < oo for n < oo. Hence, Qfjn is finite because the sum in 1)38]) 
has at most In — 1 non-zero terms. 

Next, we consider the probability distribution, P m (z'\z), which is the probability 
that a site having z particles which is charged m times, is left with z' particles after 
s = z' — z + m topplings. In the stationary state, this is related to P(s; 1, m) by 

n-1 

P(s; l,m) = } j p z Pm{m + z - s\z). (39) 

2=0 

Since the values of z a site passes through as it topples is a Markov chain, we may 
write down the Chapman-Kolmogorov equation for P m (z'\z), 

n-1 

P m+l (z'\z) = Pi(z'\z")P m (z"\z). (40) 

2" = 

The probabilities Pi(z'\z") are related to the single site operator, Gi by Pi(z'\z") = 
(e z i\iGi\e z ")i and so we introduce the matrix P m , such that P m (z'\z") = 
(e 2 '|iP m |e 2 »)i and 1)10]) becomes 

n—1 

(ez'|iP m +i|e2)i = y](e z /|iG 1 |e z »)i(ey/| 1 P w |e z ) 1 

z"=0 

= (e z/ \ 1 G 1 P m \e z ) 1 . (41) 

Thus P m +i = GiP m and, since Gi is a regular Markov matrix, there will be a unique 
stationary distribution P^, satisfying 

Poo = GxPoo. (42) 

The only non-trivial solution to ()42)) is 

Poo = (10)1,10)!,..., |0)x) (43) 
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where we have noted that the stationary state, |0)i, is unique and conservation of 
probability requires Y^~=o( e A~^<x>\ e z)~L = 1 This leads us to 

lim P m (z'\z) — p z > (44a) 

m^oo 

lim P(s; l,m) = y^ p z p {m - s)+z (446) 

z 

where we have noted that (e 2 /|0)i = p z i and the sum is over < z < n — 1 satisfying 
0<m-s + z<n-l. Finally, ()446|) leads to 

n-1 

lim Q?i = s2 Ep^ P«- ( 45 ) 

m— >oo ' ' • ' • 

s=— n+1 « 

where the sum is over all < z < n — 1 such that 0<z + s<n — 1. This limit 
exists, is non-zero and is finite for n < oo. 

Finally, we show that Q\ m > for all m. Consider again the probability 
distribution P m = [Gi] m and the equation for the width, 

oo n— 1 

<5il = S( s -™) 2 5Z P* Pm(m + z-s\z) (46) 

S = 2=0 

~ (21 

For ^* = 0, for some m*, we require 

Pz(e z >\iP m *\e z ) 1 = p z 5 z>z > (47) 

which follows from normalisation of P m and the fact that ()46|) has no negative terms. 
This implies that, for all < z < n — 1 for which p z > 0, 

[GiT* |c.)i = |c»>i (48) 

where iV > is an integer. However, since Gi is regular, there exists an integer iV* 
such that there is only one vector, |0)i, satisfying 

[Gi]>)i = |0)i (49) 

for any N > 0. If there are more than one values of z for which p z > this contradicts 
pSjl . and so Qi m is never zero. If, however, we have a single value, z* such that 
= $x,z*i then 1 0) i = \e z *) and PHj) does not lead to a contradiction. However, in 
this case the dynamics are trivial as the steady state has all sites with exactly z* 
particles and any particle added to the system will pass through immediately with 
exactly L topplings. 



5. Mapping To Area Under Random Walker 

We now come to the main result we need in order to determine the avalanche 
statistics for the directed sandpile, which is that it may be mapped exactly onto 
a random walker on [0, oo) with an absorbing boundary at the origin. After adding 
a particle at the beginning of a time step, site i — 1 will topple si > times with 
probability P(si; 1, 1). These si particles are redistributed to site i = 2, which will 
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topple S2 > times with probability P(s2] 1, si). The probability of site 2 toppling 
S2 times, independent of si, which we denote 02 (£2), is 

00 

2 ( S2 ) = ^P( S2 ;l, Sl ) (50) 

81 = 1 

which follows from ([36)1 . Defining as the probability that site i topples 2 times 
independent of previous topplings, we have 

00 

<j) i+ i(x) = ^2<j)i(y)P(x; l,y) for i = 1, . . . , L - 1. (51) 
j/=i 

This is a random walker on the interval [0, 00) with the probability of hopping from 
y to x equal to P(x; l,y). There is an absorbing boundary at x = since any non- 
toppling site stops the avalanche. If we denote the trajectory x(i), i = . . . L, then 
the avalanche size is 

L 

s = X>W (52) 

i=l 

with x(0) = 1, which is the area under the trajectory x(i). 

Note that the random walker described by (fBTj) has jumps which are correlated 
since the probability of hopping from y to x depends explicitly on y and x, and 
not simply the difference x — y. This means we must be careful if we wish to use 
the results for the uncorrelated random walker, or its continuum limit. However, 
in Ref. [TT] . the author remarks that for martingales with a fixed maximum jump 
size exhibiting stationarity and ergodicity, there is a quantity, sf = E Ylli=i °n sucn 
that 

\imP[x(i)/ Si < x] = (2tt)" 1 / e"^ 2 dy (53) 
t-* 00 J-00 

where a n 2 is the variance of the nth step in the process. 

To apply this result, we extend the random walker described by P(x; l,y) to 

the full space, (—00, 00) for we may add in the effect of the boundaries later by use 

of mirror charges 15]. As we have assumed the existence of a unique stationary 

state and have proven that Qx m = m and < Qf'm < (n — l) 2 , all that is left to 

prove is ergodicity. This is equivalent to showing that the set of recurrent states of 

the random walker are irreducible, that is, the probability of reaching any recurrent 

state i from any other recurrent state j is non-zero. Two states, % and j which 

have this property are said to intercommunicate, denoted i «-> j. We consider the 

fact that Gi is assumed to be regular, in which case there exists an N such that 

(e 2 |[Gi] m |e 2 /) > for all z,z' G [0,n — 1] and m > N. Hence, all states i,j > N 

intercommunicate since P(m ± 1; 1, m) > for all m > N. We also note that 

<-> N and 1 <->• iV which follow respectively because the avalanche should always 

be able to finish in an infinite system and arbitrarily large avalanches can be initiated 

from a single added particle. When we consider states k < N, we note that there 
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can only be a finite number of these which do not intercommunicate with state 
1. Since there is a unique stationary state which includes all states i > N, these 
non-intercommunicating states must be transient and ergodicity follows. Hence, we 
have now proven that for a toppling rule obeying the restrictions (i)-(iv) with a 
unique stationary state, (i.e Gi is regular), the distribution of the random walker on 
(—00, 00) will approach the normal distribution given by (|5Hjl . This means that, for 
long times, such a random walker with dependent jump sizes will have the statistics 
of ordinary diffusion with diffusion constant 2D = s n . Hence, by adding mirror 
charges to remove paths that cross x = 0, we are able to calculate the large L 
statistics of avalanches directly from the area under the Brownian curve, which is 
our justification for calculating moments in the continuum limit in the next section. 
Of course, we could have simply gone ahead and carried out the calculations in 
the continuum without the above analysis and demonstrated that they correctly 
modelled the numerics. However, had we done so we would not have had a precise 
idea of how trustworthy these calculations were and where we expect them to break 
down. 

6. Moments of the area under the Brownian curve 

Having proven the correspondence between avalanches and a random walk of 
independent identically distributed step sizes, we proceed to calculate the moment 
generating function for the area under the Brownian curve. The authors are aware 
of only one study which investigates the finite-size effects due to stopping the curve 
after some time, which corresponds to the finite size of the sandpile [2j and since 
our analysis goes further than that in Ref. [2], we present it here in some detail. 
The following calculation will be carried out using notation and language suitable 
for the random walker description of the problem. Hence, the Brownian curve will 
be described by a trajectory x{t) where x is interpreted as "space" and t is "time" 
with the diffusion constant D having units Length 2 Time _1 . We do this because the 
path integral approach we are about to employ is more intuitive in this language. | 
The connection to the sandpile is made by noting that the number of topplings of 
site i is equal to x(t — i) and the system size, L, is equal to time at which we stop 
the curve x(t). 

We begin with the generating function 

POO 

(e- XA )= e- XA P(A;x,L)dA (54) 
Jo 

where P(A; x, L) is the probability that a random walker starting at x has the area 
under its trajectory equal to A after time L. If we denote the trajectory of the 

I We should stress that in this picture we consider the Brownian curve x(t) as existing on the 
entire interval [0,oo) and we measure the area up to the point L, A = J Q x(t)dt. Hence, what is a 
boundary in the sandpile picture (the open boundary at site i = L) is not considered a boundary 
in the Brownian curve picture. 
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x(t) 




L 



Figure 1. The area under the Brownian curve with an absorbing boundary. The 
statistics consist of contributions due to curves x(t) which are nonzero at t = L 
(solid line), as well as those which cross the boundary at some time t < L (dashed 
line). 



walker x(t), then the curves contributing to P(A;x,L) are all those which satisfy 
x(t)dt = A. (55) 



Note that we have an absorbing boundary at x = 0, such that if x(t') = for any t 
then x(t > t') = 0. Hence, there are two contributions to P(A;x, L): That due to 
trajectories which do not cross the absorbing boundary, x(L) > and those which 
cross the absorbing boundary at some time t < L, see figure |H1 
We shall treat these separately, writing 

roo /*oo /»oo i*L 

( e - XA )= dA dy e~ XA ^(A,y,L;x) + dA dt e" AA $(A, t; x) 
Jo Jo Jo Jo 

= h + h (56) 

where ^(A, y, L; x) is the probability that a trajectory beginning at x(0) = x passes 
through x(L) = y with area A and §(A,t;x) is the probability that a trajectory 
beginning at x(0) = x first touches the absorbing boundary at time t, with area A. 
Using standard path integral methods, we may write down 

V(A, y, T; x) = lim / &x{t)8 ( [ x(t)dt - a] exp (- [ [Dx 2 + gd(x))dt) (57) 

9-°°Jx(0)=x \Jo J \ Jo J 

x(L)=y 

where x = dx(t)/dt. Taking the integral over A we find that the first term on the 
right hand side of (JHEJ) is 

h = lim [ dy [ @x{t) exp (- [ \Dx 2 + Xx + g5(x)] dt] .(58) 

9^ Jo Jx{0)=x \ Jo J 

x(L)=y 

Following the lines of Ref. [12] . we note that this is simply the path integral for a 
Brownian particle with a linear potential for x G (0, oo) and an infinite potential at 
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x = 0. Hence, we write this term as 



Ji = lim 

s^ 00 Jo 



-HL\ 



x)dy 



(59) 



H\ 



where H = —D-jj^ + \x + g5(x). The resulting equation of motion, J^|0) 
is easily solved using Airy functions which can be used to form an orthonormal basis 
on [0,oo) [IH], 



Aife^x + Xj) f x °°M(z)dz x 



A 2/ 3jD l/3 L 



where Xj are the zeros of the Airy function, x\ = —2.338 . . ., x 2 = —4.087. . . etc. 
In a similar way, for the second term on the right hand side of (JHEJ), we have 



(60) 



d 

$(A,t;x) = D —V(A,y;x,t) 
dy 



(61) 



since this is the current of diffusing particles with area under the curve equal to A, 
leaving the system at time t. Hence 



L 3 
o dy 



dt 



Ai ((b) 



AN 1/3 



X + X 



Ai'lx-i 



(62) 



3=1 

In order to proceed, we use the fact that the leading order L dependence for 
each moment come from terms linear in x. In Ref. [2] it was shown that if the 
moment generating function is written, 

gn 



(A n ) 



-l) n 



A=0 



= n\ <f> n (x,L), 
then <f> n (x,L) may be determined recursively 



<f>n{x,L) 



roc pL 

I dx' / dt x'(j) n -i(x', t)G(x, L — t; x'), 
Jo Jo 



(63) 



(64) 



where G(x, t; x') is the propagator for the diffusion equation with appropriate 
boundaries, 



G(x, t; x') = 

If we define the "current" 

d 



(x — x ) 

e~ iot 



(x + x') 2 

e iot 



V^irDt 



3n(L) = —(f) n (x,L) 



(65) 



(66) 



x=0 



then if j n (L) is non-zero, (j) n (x, L) is proportional to x to lowest order. In this case 



Jn{L) 



dx' 



dt 



X 12 4>n-l{x',t) 

4^[D{L-t)fl* 



g AD(L-t) 



(67) 
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B and so j n {L) > for n > since the integrand is always positive definite. Hence, 
all moments are proportional to x to lowest order. The fact that the terms linear 
in x will also be the highest order in L follows from dimensional analysis. If we 
write down the expansion of a moment in powers of x then each term must have the 
same dimension. By considering the dimensions of the available quantities, such an 
expansion must take the form 

oo 

(A n ) =xC„D("- 1 )/ 2 L( 3n - 1 )/ 2 + ^x fc C'„, fcJ D("- 1 - 2fc )/ 2 L( 3 "- 1 - 2fc )/ 2 (68) 

k=2 

where C n> k are simply more coefficients with no x, D or L dependence. Hence, the 
term of lowest order in x will have the highest order L dependence. 
Taylor expanding 1% and I2 to first order about x — 0, 

h « xX / V e^ x2/3Dl/ ' H dt = J 2 . (696) 
Jo J=0 

Note, however, that this approximation is not valid for the zeroth moment, (A ) = 1 
since it is not proportional to x. J\ and J2 are now in similar forms to equations 
appearing in Ref. [T2j. They calculate the quantity 

~ ./LT Ai( * )ck 
Ma*; 



P(X,L) = v / 7t2~ 1/6 (AL 3/2 ) 1/3 Y Jx \ , , e ^ A2/32 ' 1/3L 



n=0 



n\ 



where the a n have been calculated in Ref. [T7] . We simply quote the first few values, 

3/7 59 465 pi: 5345 

°° = 1 ' ai = 4V2' a2= 60' a3= 512V2' a4 = 3696- 
Apart from a few multiplicative prefactors, (|70p differs from ()69aj) only by the fact 
that the former uses D = 1/2. We therefore have to reinsert the diffusion constant, 
D, which they assumed equal to 1/2, but this is easily done by considering the 
dimensions of the results. We note that J\ is a dimensionless function, and so 

y^'r^HKL) (72) 

where 7 is chosen such that J\ is dimensionless. It is then easy to show that 
7 = 1/6 — n/2 and 

J 1 = x y iZ^! CnjD («-l)/2 L (3n-l)/2 (73) 





where 

2 n / 2 a r 



(74) 
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We may carry out an identical procedure for J 2- The equivalent quantity in 
Ref. (H is 



P(X,L) = V2^(XL 3/2 )J2 e 



oo 

, A 2/3 2 -l/3 L 



3=0 



£ ^~b n L^ (75) 



n 

n=0 



where, again, b n have been calculated in Ref. [T7j, the first few values being 

1 F , 5 15 /jt , 221 

On = 1, 0i — —\—, 02 = , Os = 1 / — , 64 = . (76) 

° ' 2\/2' 12' 3 64 V 2' 1008 V ; 

Following the same steps as above we find 

j 2 = x y (z^lV 2 -^/^- 1 )/ 2 [ L &»-Wdt 

^ n 20F Jo 

n=0 V JU 

Y iz^l dnD (n-D/2 L (3n-l)/2 ^ 



X 

n=0 



where 



dn = Tl \ \ (78) 
(3n - 1JV7T 

Hence we have 

(e~ XA ) = 1 + x Y C„^4^^ (n_1)/2 ^ (3 "" 1)/2 + ^V) (79) 
^— ' n! 

n=l 

where C n = c n + and the first few values are 

Ci = 1, Co = Ca = — , Ca = ;=. (80) 

3 8 ' 693^ v ; 

The first two values are in perfect agreement with those derived in Ref. |2j, and the 
authors are unaware of any previous calculations of C n for n > 2. Thus we may 
immediately identify the exponents r = 4/3 and A = 3/2 and the amplitudes allow 
us to compute universal amplitude ratios, which we will use later to compare the 
numerics against theory. 



6.1. Crossover from Branching Process 

The convergence of <pi{x) to the normal distribution occurs only as i — > oo, and 
hence the results above are only valid for L — > oo. In using the Brownian curve 
instead of the exact curve described by P(s; l,m) we have taken a hydrodynamic 
limit and therefore thrown away any information about the statistics of the process 
for small L. It is natural, therefore, to ask how we expect the results to differ in this 
regime. We propose the existence of an n dependent crossover length, £ n , such that 
the above scaling analysis is valid for L ^> £ n . We argue that for smaller systems, 
1 <C L <C £ n , we expect to see scaling corresponding to the branching process. 



Id Directed Sandpile Models and the Area under a Brownian Curve 



17 



Consider adding a particle to the first site. If the probability that the site has 
z particles, p z , has support for all z G [0, n — 1], then for n 1 it is likely that 
0<z<Cn-l. In this regime we may assume that the number of times the site 
topples due to this added particle, which we denote si, is largely independent of 
z. Site 2 therefore receives s\ particles, each of which may cause it to topple s 3 2 
times, i — 1 . . . s% with the total number of topplings of site 2, S2 = s 3 2 . While z 
remains far from and n — 1, the s J 2 will be largely uncorrelated and by continuing 
this argument to more sites, we see that while Sj remain small each site will topple 
nearly independently. However, as we continue through the system to higher i the 
Si will begin to see large fluctuations and the avalanches will become correlated, 
assuming the scaling of the previous section. Hence, we argue that for systems 
with 1 C L < £ n , the avalanches will resemble those of the uncorrelated branching 
process with exponents r = 3/2 and A = 2 [IF] . For larger systems, L ^> £ n , 
temporal correlations emerge in the avalanches and r = 4/3, A = 3/2. 

The fact that the above argument relies on realisations where p z has support 
for a large range of z indicates that the crossover length £ n depends on the details of 
the toppling rules and as such cannot be thought to have any "universal" qualities. 
Indeed, we have not specified how the toppling rules in a realisation should be altered 
as n is increased, and so it is impossible to say anything a priori about the behaviour 
of £ n - 

7. Numerics 

We now support our claims with numerics by demonstrating that the correct scaling 
(with crossovers - see previous section) occurs for a particular realisation of this 
directed sandpile model. In order to study the scaling we choose a realisation such 
that it is clear how to generalise to higher n. The only remaining difficulty is to 
find the correct variance to put into the equations when we come to compare with 

~ (2) 

numerics. In all that follows, we use 2D = Q\ t00 as we find that it fits the data very 
well. 

We compare the scaling predicted above with numerics from a realisation with 
the following toppling rules: A site i, 1 < Z{ < n — 1 which receives a particle will 
topple 1,2 or 3 times with probability 1/8 or will not topple with probability 5/8. A 
site with z^ = will topple once with probability 3/8, a site with z^ = 1 will topple 
once with probability 2/8 and twice with probability 1/8 and a site with Zi = n — 1 
will topple once with probability 6/8, and 2 or 3 times each with probability 1/8. A 
site with zi = n — 1 has to topple at least once in accordance with restriction (i). 

We expect (s 2 ) to scale with the system size 

(s 2 ) l ~{l l«L«( n , L 5 / 2 L>£„, (81) 

where £ n is a correlation length with some (as yet unknown) n-dependence. These 
results have been confirmed and are shown in figure El 
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Figure 2. Numerical results for n = 4, 8, 16, 64. The errors for both graphs were 
calculated using Efron's Jackknife |T2|. (a) The rescaled second moment (s 2 )l 
vs. system size. For large systems this is a constant for all values of n. (b) The 
moment ratio gz{L). For large L this approaches the constant value g$ ~ 1.29 for 
all values of n. 



We also analyse the moment ratios defined by 

a (L) - < a *M*>E" 2 (82) 

It is a straightforward calculation to show that, for an avalanche probability given 
by (JIJ gk{L) approach universal values for L — > 00. These values are simply ratios 
of the amplitudes C n calculated in |Hl 

g k = lim g k (L) 

L— »oo 

( (83) 



ct 1 



This agrees with the numerics, as illustrated in figure 01 for gs(L) which appears 
to converge to a universal value of 53(00) ~ 1.29, in excellent agreement with the 
theoretical prediction g% = 15 3 7r/2 13 = 1.294. . . as well as numerics for a different 
realisation published elsewhere ^U]- This supports our claim that the L — ► 00 
limits of the gk{L) are indeed universal. Note also that £ n has a notably different n 
dependence in this model than in the one presented in [TOj- In this case, £ n saturates 
to a constant value for large n, meaning that for large n the crossover occurs at the 
same value of L. This is because the support of p z is finite for n — » 00. 



8. Conclusion 



We have found the stationary state avalanche-size distribution for a general n-state 
directed sandpile model. The avalanches can be mapped onto a random walk of 
dependent random variables and, using an applicable central limit theorem, we have 
shown that under a broad set of conditions the moments scale with r = 4/3 and 
A = 3/2. We also note that this value of r agrees precisely with that obtained in 
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Figure 3. Numerical results for n = 4, 8, 16, 64. The errors for both graphs were 
calculated using Efron's Jackknife ^H] and are approximately the same size as 
the symbols, (a) The moment ratio gs(L) (a) and 54 (L) (b). For large L these 
approach the constant values g% fts 1.29 and 34 s» 1.9 respectively for all values of 
n. The dashed lines indicate the exact values <?3 = 1.2942 . . . and 54 = 1.8975 . . ., 
in excellent agreement with the numerics. 

[T3] . which calculates the probability distribution in the infinite system size limit. 
We have also calculated the moment generating function for the area under a random 
walker with an absorbing boundary, and found a relation for the moment amplitudes 
in terms of those already known for other Brownian processes. 
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Appendix A. Explicit Example: n — 2. 

In this section we calculate the steady state properties of an n = 2 model, which is a 
generalisation of the model studied in Ref. |2j, and compare predictions to numerical 
simulation. For n = 2, the most general model we can write down, which obeys the 
rules (i)-(iv), is 



where < a < 1 is the probability that a site with z = does not topple on receiving 
a particle and < (3 < 1 is the probability that a site with z = 1 topples twice on 
receiving a particle. Hence, the single site toppling matrix is 



We now proceed to calculate the steady state properties of this model, following 
the prescription given in El Gi(x) has eigenvalues A = x, \i — x(l — a — (3) and 
eigenvectors 




(1.1a) 



(1.16) 



(1.1c) 




(1.2) 




(1.3a) 



(1.3d) 



(1.36) 



(1.3c) 



Hence, the eigenvector for the stationary state is 




(1.4) 



valid for ^ 1. 

From these results it follows immediately 



~ (2) a 

Qi.L =PoPi+PiPo = 2 — 

a - 

and hence, from (f7U|) and using 2D = Q 




(1.5) 
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Figure Al. Numerical results for Q\ ^ for a — 1 as a function of p with data 
(squares) compared against the values predicted by l|1.5|l (solid line). The data 
was obtained by measuring Q\ ' for large m and estimating the asymptotic value. 
Comparison is made across the whole range of (3 and the inset shows data in the 
vicinity /3 — > 1. Note that the agreement is excellent right up to = 1. Typical 
error bars for the numerical data are the size of the squares. 



in perfect agreement with numerics, see figure lA"T1 and figure [XU However, it should 
be noted that for a and f3 both approaching 1 the random walker it describes will 
spend more and more time on either only odd or only even sites. Hence, it will take 
longer times (larger system sizes) for the statistics to reach the asymptotic values 
and so we expect very strong corrections to scaling for a, [3 — > 1. When a — (3 — 1, 
we no longer have a unique stationary state and so scaling is not observed. 



Appendix B. Calculating the Amplitudes C n 

The amplitudes a n and b n appearing in El can be calculated using the methods 
outlined in Refs. O E] . For a n we define 

, ^ „r(l/2)n! , x 

On = (^Y^pyRn (2-1) 

where R n are constructed through the following recursion relations 

n 

R n = (3 n -Y,lj R n~3 (2-2) 

3=1 

p n = 7n + |(2n-l)/3 n _i (2.3) 



4 

r(n + l/2) (36) n n!' 



7,. ^ ^i±J#^. (2-4) 



Similarly, for the b Ti 



b n =A(V2)~ n ^^K n (2.5) 



2 / 
n-1 



K n = -^—K^ + K 3 K n-r (2.6) 



3=1 
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(a) 




(b) 




Figure A2. Numerical results for the rescaled second moment (s 2 )l/L 5 / 2 for 
a = 1 (a) as a function of j3 for L = 256 (*), 1024 (o), 4096 (o), 16385 (A), 65536 
(□) and 131072 (x). (b) Rescaled second moment (s 2 )l/L 5 / 2 for /3 = 0.95 vs 
inverse system size. The dashed line is the theoretical value. The measurements 
appear to converge towards the theoretical line large L, supporting our claim that 
the deviation is a finite size effect. 



Putting these together and rearranging slightly, we find that the amplitudes C n are 
given by 

Cn = f ^iy (Rn + 2K n ) (2.7) 

We tabulate the first 10 values of C n , along with the universal amplitude ratios 
9n = CJC^' 1 in table Bl. 
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Table Bl. Tabulated values of C n and g„ 



On 



1 1 

2 
3 



4 
■5 
6 
7 
8 
9 

10 



:v> 



:S:S75 



7T 



15V¥ 

15 

8 

4096 
693^ 

2875 
448 

1219336 
51051^7? 

745039 
24576 

25796624240 
200783583^ 

214422265 
1171456 

15033906553126 

17468171721^ 3793868231748622483456 ' 



8192 

47625 , 
78848 ' 



145546875 2 
469762048 



71374471168 



2828819953125 3 
8796093022208 77 



10202766423046875 3 
15969609677012992 



549540812759765625 4 
1288029493427961856 
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